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Abstract 

A previous study of diatomic molecules revealed that variational second-order density matrix 
theory has serious problems in the dissociation limit when the N-representability is imposed at the 
level of the usual two-index (P, Q, G) or even three-index (Ti, T2) conditions [H. van Aggelen 
et al., Phys. Chem. Chem. Phys. 11, 5558 (2009)]. Heteronuclear molecules tend to dissociate 
into fractionally charged atoms. In this paper we introduce a general class of A^-representability 
conditions, called subsystem constraints, and show that they cure the dissociation problem at little 
additional computational cost. As a numerical example the singlet potential energy surface of 
BeB"*" is studied. The extention to polyatomic molecules, where more subsystem choices can be 
identified, is also discussed. 
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I. INTRODUCTION 



In recent years much attention has been devoted to the direct determination of the second- 
order density matrix (2DM) through variational optimization. As first shown by Husimii, 
the energy of a system interacting with at most two-particle interactions is fully deter- 
mined by the 2DM. Some fifteen years later Lowdin^ independently derived similar results 
and suggested determining the 2DM directly in a variational approach. The first practical 
calculation was done by Mayer- who tried to compute the energy of an electron gas by a 
variational optimization of the 2DM. The energies obtained however, were much too low, 
and inconsistent with existing results. Tredgold^ realized that the problem arises because 
the set of density matrices over which the optimization is carried out is too large. Although 
in these first attempts some obvious constraints were included, better constraints are needed 
in order to make sure that the 2DM can be derived from a physical wavefunction. This prob- 
lem was termed the A^-representability problem by Coleman in his seminal review paper-, 
in which he solved the ensemble A^-representability problem for the first-order density ma- 
trix (IDM) and derived some bounds for the eigenvalues of the 2DM. Garrod and Percus^ 
subsequently derived the much stronger positivity conditions Q and G. Because of the 
computational complexity and some dissapointing results on nuclei^, not much progress was 
made the next twenty-five years. Interest renewed in the direct variational determination 
of the 2DM after Nakata et al.- and then Mazziotti^ used a semidefinite program algorithm 
(SDP) to study a number of small atoms and molecules and got reasonably accurate re- 
sults. These results sparked of a lot of developments. New A^-representability conditions 
where introduced, e.g. the three-index T conditions, as set forth by Zhao et al.— , which 
led to mHartree accuracy^"— for molecules near equilibrium geometries, and generaliza- 
tions thereof — . Algorithmic breakthroughs were realized with the implementation of 
a scaling SDP algorithn>22"^ and the development of an active-space variational 2DM 
method^"—. A drastic failure of the standard A^-representability constraints {PQGT) was 
shown to occur in the dissociation limit by van Aggelen et al.— using a recently developed 
semidefinite programming code^^. In this article we propose new strict constraints, which 
we call subsystem constraints, that fix the inaccuracies in the dissociation limit. Sec. [TTl 
contains the theoretical derivation of the subsystem constraints in a general framework. It 
is shown that identifying a subspace of the complete single-particle space leads to upper 



2 



bounds for the energy of the total system, that must be obeyed by any A^-representable 
2DM. As a simple illustration, the technique is applied in Sec. IIIII to the dissociation of 
BeB"*" in a small (Dunning-Hay) basis set. Sec. HV] contains a summary and discussion. A 
systematic and thorough study of the diatomic potential energy surfaces for the 14-electron 
series is the subject of a separate publication^^. 



II. THEORY 



A. Integer-A^ ensemble representability 

The second-order density matrix (2DM) corresponding to an A^-fermion wavefunction 
1^-^) is defined as 

Second-quantized notation is used where aj^, [aa) creates (annihilates) a fermion in the 
single-particle (sp) state a. The sp basis is assumed to be orthonormal throughout the 
article. Eq. ([1]) is easily generalized to the 2DM corresponding to an ensemble of A^-fermion 
wavefunctions; conversely, a matrix is called integer- A^ ensemble representable if 

r^ft,5 = J]x,(vl;f |alata,a,|vl/f) (2) 

i 

for some set of A^-fermion wave functions and positive weights Xi obeying = 1- 

We consider a system governed by a Hamiltonian containing a one-body part t and a 
two-body interaction V, 

H = ^ta^al^a^ + ^^Vap.^sala^pasa^, (3) 

where Vap-^^ys represents the antisymmetrized matrix elements of the interaction. The exact 
ground-state energy (assumed to be nondegenerate) is determined by finding the 2DM 
that minimizes the energy functional 

= mill [Tr(tp^) + Tt{VT^)] , (4) 
is integer- A^ ensemble representable 
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where is subject to the constraints hsted after the {-symbol. Note that the Tr(ace) 
operation in the 2DM space is restricted to antisymmetric two-index combinations and that 
the first-order density matrix is defined through the partial trace = Tq,^.^^. 

A lower bound to the exact energy is obtained by the variational determination of the 
2DM subject to a selected class of necessary A^-represent ability conditions. This reduces to 
finding the 2DM that minimizes the energy functional 

E§nP = min [Tr(tp^ + Tt{VT^)] (5) 

< TrT^ = |iV(iV- 1) 
Ci{T^,N) > 

The Ci{r, N) are matrix functionals of the 2DM which are required to be positive semidefi- 
nite; they reflect a choice of necessary conditions for integer- ensemble representability, e.g. 
the two-index P, Q and G condition or the three-index Ti and T2 conditions. As indicated 
in the notation, the matrix functionals also depend on the particle number. 

B. Fractional-A ensemble representability 

The exact solution for a fractional electron number A" has only one sensible definition, 
based on considering ensembles containing wave functions with various electron numbers, 
and the resulting energy then has the well-known piecewise linear behavior between integer 
valued'—. 

Reformulated in terms of density matrices, one defines a 2DM F to be A^-representable if 

F = ^x^F^ (6) 

TV 

for a set of integer- A^ ensemble representable F^ and a set of positive weights obeying 
Y^j^Xn = I and Xljv = N. 

Obviously, the IDM p corresponding to the same ensemble cannot be obtained directly 
from F by a partial-trace operation without knowledge of the ensemble weights. It is therefore 
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more natural to consider the pair {p, T) as being fractional- ensemble representable if 



ViV : = ^ 
^ VA^ : is integer- ensemble representable 

The exact solution is then simply generated by the minimization problem 



- min [Tr(tp) + Ti{VT)] 



(p, r) is fractional- ensemble representable 



(7) 



(8) 



It is now clear how the variational problem for a selected choice of N-representability 
conditions [corresponding to Eq. (^] should be phrased, when it is generalized to a fractional 
electron number: one should minimize 



E§^P = min [Tr(tp) + Tr(V^r)] 



(9) 



p = T.N^Np^;'^ = Etv^w^^ 

VAT : p^ = 

VA^:Trr^ = |A^(A^-l) 
^ VA^ : Ci{T^,N) > 

where both the weights xn and the 2DM's can be varied. For any choice of weights 
the energy is minimal when corresponds to the SDP solution for integer A^, so Eq. (Q 
can be reformulated as 



EsDP = mm 



'SDP 



Y^x^Ei 

N 

{ a:Ar > 0; Y^j^xn = 1; Y.N = N 



(10) 



This leads naturally to a piecewise linear solution^^ which, for a convex set Egj^p, is given 
by 

EinP = {N- IntmE'sSP"-' + (Int(iV) + 1 - N)E';^f^ (11) 

where the function Int(a;) returns the nearest integer number smaller than x. So in exactly 
the same way as for the exact solution, nothing new emerges from introducing fractional 



electron numbers. There is, however, one interesting observation, that is crucially important 
for the subsequent derivation of subsystem constraints: for any Hamiltonian [t, V) the energy 
evaluated with a fractional- iV ensemble representable density (p, F) obeys the inequalities 

E{p,T) > > E§^p. (12) 

The first inequality follows from the definition (|9]) of Eq as the minimum over ensembles, 
the second one from the fact that for all integer N, Eg^p is a lower bound for E^ . 

Some electronic structure methods (like Hartree-Fock or density functional theory) can be 
easily generalized to accomodate a fractional number of electrons. In that case the behavior 
of the energy in between integer values usually disagrees with the piecewise linear result 
from an exact calculation^"—. 

Something similar is also present in the SDP technique. A naive extension of this frame- 
work to fractional electron number would be to treat the parametric dependence on in 
Eq. ([5]) as a continuous variable: this is immediately understood to be unphysical, as no 
reference is made to the ensemble interpretation. However, as first noted by van Aggelen et. 
al.— , SDP computations on the union of isolated subsystems become equivalent to allow- 
ing the number of electrons on each subsystem to be a continuous variable, again without 
reference to the ensemble. Application of the SDP to the dissociation limit of diatomic 
molecules and ions results in dissociated atoms with fractional occupancies in the case of a 
heteronuclear diatomic, which is chemically clearly unacceptable. In order to prevent this, 
new A^-representability conditions are needed. 

C. Subsystem constraints 

As before, Greek indices a,/3, .. are used for the full set of single-particle states. We 
now introduce an (arbitrary) subset; Roman indices a,b,.. are used if we want to restrict 
the orbitals to members of this subset. In a polyatomic molecule, e.g., we may consider an 
orthonormal basis of the subspace generated by the basisfunctions centered on a particular 
atom. This selection of a subset is in fact equivalent to choosing a subspace of the complete 
single-particle space, since the SDP setup with the standard two-index or three-index con- 
ditions (of the P, Q, G, Ti, T2 type) is invariant under a unitary transformation in sp space. 
Note that this invariance is broken when in a particular sp basis only diagonal constraints 
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are kepUii^a-is, 

The derivation of the subsystem constraints starts by noting that, if is integer- 
ensemble represent able, then the pair (p*^"^, F*^"^), 



sub 



Pac 



—^T.^a,;c, = pl (14) 



is fractional- ensemble representable in the Fock space generated by the subset, with 
N = YlaPaa"- '^'^ prove this we start with the integer- ensemble representability of F^: 

i 

We can expand each {"^f) in Slater determinants, classified according to the number of 
subsystem orbitals they contain: 

l^f)= E C.3s,..~Ms^-3) (16) 

in which < j < A^, sj represents a set of j subsystem orbitals, and sn-j a set of 
N — j orbitals not in the subsystem. Using the fact that the string of subsystem-type 
creation/annihilation operators in Eq. ( |T5l) does not change the number of subsystem or- 
bitals, that it leaves the non-subsystem part of the Slater determinant unchanged, and using 
orthonormality of the sn^j states, we see that 

i i jsN-j 

where 

i*L._; = E^^^.^-->^) (18) 

is a state with j particles in the Fock space generated by the subsystem orbitals. These 
states are obviously not normalized, their norm is given by 



C^h^jn.J = E \C^3S^^N-f = ■ (19) 



If we replace them by normalized states 

1^^- ) = \wt ) , (20) 



it follows that 

j;iSN-j 

where 

E ^^^L^. = 1 ' (22) 

i;isjv-j 

because of the normalization of the original A^-particle states. In an analogous way one 
shows that the first-order density matrix p^^^ can be written as 

p:t= E ^^<.-M.J^^M^J■ (23) 

This proves that T^^^ and p^^^ can be derived from the same ensemble of wave functions 
containing only orbitals in the subsystem. This ensemble has a fractional number of particles 
(in the subsystem space) given by = 'Zj.isj,., j^Msj,., = J2aPTa- 

Based on Eq. (fT2|) . the following necessary condition then holds: If is integer- 
ensemble represent able, then (p*^^*^, F'^^*') should obey 

Tr (t^'^^p^^*^) + Tr (y^^bpsub^j > ^n^^ ^24) 

with N = XlaPa^' ^'^y Hamiltonian (t^^^^, 1/*^"^) defined in the subspace. These subsys- 
tem constraints can be quite powerful; they are grossly violated in Coulombic dissociation 
problems as documented in^^. Note that for consistency it is preferable to use the SDP lower 
bound for the subsystem energy rather than the exact one, even if this should be available. 

The subsystem constraints in Eq. ( l24l) were derived for a particular basis choice of the 
subspace, i.e. a subset of the underlying orthonormal basis \a) of the total single-particle 
space. However, the subsystem constraints only depend on the subspace itself, which is most 
easily seen by extending the t^^^ and V^^^ operators (defined on the subspace) to the total 
single-particle space using projection operators. One can rewrite, e.g., 

Tr (t-V^^') = E = E /^"7^?' (25) 

ac 07 

where 

C7 = (a|Pf"^P|7). (26) 
Here f""^ and P are first-quantized operators with 

P = J2\a){a\ (27) 

a 

8 



the orthogonal projector onto the subspace. In the same way one can rewrite the two-body 
operator defined on the subspace as 



D. Implementation for diatomics 

For a diatomic molecule the molecular SDP solution tends to localize a fractional number 
of electrons on a well-separated atom, whenever this situation is energetically favorable in the 
continuous- sense mentioned at the end of Sec. Ill B[ As a result, the energy of the atomic 
subsytem drops below the true (ensemble-based) energy and a violation of the inequality 
(12^ occurs. In this case the choice of the subsystems is therefore obvious: in the dissociation 
limit they should coincide with the individual atoms. While a general basis-set independent 
formulation is possible, in practice the calculations are performed using atom-centered basis 
functions for which this requirement is automatically satisfied. 

The procedure for applying the atom-A subsystem constraint in a diatomic AB can then 
be summarized as follows: 

1. Solve the atomic SDP problem for a central charge Za at various electron numbers, 
using the sp orbitals centered on A. In practice, only electron numbers near atomic 
neutrality are important. This generates the atomic energy E'^ as a function of frac- 
tional electron number N (i.e. a piecewise linear curve). 

2. For each internuclear distance Rab-, calculate the transformation matrix between the 
(nonorthogonal) atomic basis of the A-centered orbitals and the orthonormal 
molecular basis |a) that is used in the SDP program. 



The U"^ matrix is easily constructed with standard quantities in molecular modelling 




(28) 



by defining 




(29) 



(30) 



packages. 




(31) 
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with C the expansion coefficients of the \a) molecular basis in terms of all the 
nonorthogonal orbitals centered on the various atoms, 



) = Y.^c.;M. (32) 



and with Si^-j^ = («a|jd) the overlap matrix for the atom-centered basis functions. 
The orthogonal projector onto the subspace spanned by the A-centered orbitals, when 
expressed in terms of the nonorthogonal basis set, reads^ 

P'' = J2(Sa%M{ja\ (33) 

ij 

where Sa is the block of the overlap matrix corresponding to the A-centered orbitals, 
and S^'^ is the inverse of this block. 

3. Perform the SDP program with the extra linear inequality: 

Tr {t^p^) + Tr (V^^F^) > ^^^=^^(1"^") _ (34) 

Here: 

(tX, = Y.(^A\i^\jA)Wf^Wf^ (35) 

ij 

(l^U = E('^^)^^-^>/7 (36) 

ij 

{V^)a^.,,S = Y.^iA3A\V^\kAlA)WtWfpW^^Wij (37) 

ijkl 

and tA is the kinetic energy plus attraction to nucleus A. The coefficients W are given 
by: 

Wt = Y.K{S-/),. (38) 

j 

The inequality flM|) is nothing but the application of Eq. fl2^ in the subspace defined 
by the sp orbitals centered on A and using the Hamiltonian of atom A. Obviously, 
atom B generates a similar inequality. 

The inequalities for A and B are expected to become important in the dissociation limit, as 
they prevent an artificial lowering of the energy due to fractional electron numbers on atom 
A and B. 
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III. NUMERICAL VERIFICATION 



It has been found previously that the 2DM variational optimization under P, Q and G 
constraints leads to chemically flawed dissociation curves where the atoms carry noninteger 
numbers of electrons even when separated by a large distance. This is especially true for 
molecular ions and persists even when including T constraints^^. In order to show the value 
of the subsystem constraints, we present the potential energy surface of BeB^, computed 
for a separation ranging from 1 to 9 A. The main interest here is a proof-of-principle of the 
fact that the new constraints indeed severely restrict the variational freedom in the SDP. We 
therefore opted for the fairly small Dunning-Hay basis^, making full-CI calculations still 
feasible. BeB"*" is a good example since this 8 electron system dissociates into Be and B"*". 
Application of P, Q and G for an 8 electron system is expected to yield energies that are 
signiflcantly too low compared to full-CI. However, at the dissociation limit, the energy of 
the molecule should be equivalent to that of isolated Be and B~^. As shown for both Be 
and B"*" the P, Q and G energy is very nearly equal to the full-CI energy, so application of 
the subspace constraints should result in much higher PQG energies than when not using 
the subspace constraints. 

Table [T] shows the energy of the molecule at different internuclear distances, computed 
at the full-CI level of theory as well as the variationally optimized 2DM energy with and 
without subspace constraints (see also Fig. [1]). 

Table [T] very clearly shows that the difference between the FCI and 2DM energies is sub- 
stantial when using only the P, Q and G constraints. Especially at longer separations the 
difference between full-CI and 2DM energies can amount to roughly 0.03 Hartree. The sub- 
space constraints succeed in reducing this error by approximately two orders of magnitude. 
As expected, the remaining error is very small because P, Q and G yield energies for the 
atomic 4-electron isoelectronic series that are very near to full-CI energies. The present new 
constraints are clearly very succesful. As Figure [1] shows, the constraints are active most 
for separations above 4.5 A. The nearer to complete dissociation, the more of the error is 
recovered by the subspace constraints. As shown by Van Aggelen et al.— , not only are the 
energies improved, also chemical observables and chemical concepts are substantially better 
for the 2DM obtained when including the subspace constraints. As an example, the Mul- 
liken population^ on the Be atom at 9 A is +0.38 when not using the subspace constraints. 
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TABLE I: Difference (in mHartree) between full-CI energy (FCI) and variationally optimized 2DM 
energy without (2DM) and with (2DM+) subspace constraints. 
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2DM 


2DM+ 
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2DM 


2DM+ 


1.25 


16.55 


16.55 


3.50 


13.30 


13.30 


1.50 


10.86 
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17.17 
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1.75 


9.56 


9.56 


4.00 


19.12 


19.12 


2.00 


9.94 


9.94 


4.50 


19.90 
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2.10 
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10.07 


5.00 


21.28 


14.07 


2.20 


10.07 


10.07 


5.50 


23.07 


9.35 


2.30 


10.30 


10.30 


6.00 


24.48 


5.27 


2.40 


10.59 


10.59 


6.50 


25.67 


2.44 


2.50 


10.77 


10.77 


7.00 


26.95 


1.26 


2.60 


10.91 


10.91 


7.50 


27.90 


0.66 


2.75 


11.30 


11.30 


8.00 


28.88 


0.53 


3.00 


12.00 


12.00 


8.50 


29.63 


0.38 


3.25 


12.56 


12.56 


9.00 


30.39 


0.37 



whereas inclusion of the subspace constraints yields a charge of 0.00, consistent with what 
it should be according to the full-CI data. 

The added constraints result in a much better description of molecular dissociation. Neither 
atom still suffers from fractional occupancy at the dissociation limit. Addition of each 
subsystem constraint does not slow down the SDP, as it adds a fairly simple linear inequality 
constraint. 

IV. SUMMARY AND CONCLUSIONS 

In this paper we derived necessary conditions that substantially reduce the dissociation 
problem found in diatomic molecules at large separations^. These conditions are based 
on fractional- ensemble representability. For any subset of single-particle space one can 
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FIG. 1: The singlet dissociation energy curve of BeB^ calculated in full-CI, and determined vari- 
ationally without (2DM) and with (2DM+) subspace constraints. 

associate a "subsystem" 2DM and IDM, which must be fractional- ensemble represent able. 
Any Hamiltonian defined on the subspace leads to linear inequalities for the 2DM of the full 
system. In the case of diatomic molecules we can associate the subsystems with the single 
atoms and consider the separate atomic Hamiltonians. Application to BeB"*" shows that 
these constraints are strongly violated in the dissociation limit. When they are imposed 
during the optimization, they lead to a significant improvement in the energy, and cure 
the pathological behavior reported In the dissociation limit of heteronuclear diatomic 
molecules the occupation numbers of the single atoms are now integer, as they should be. 

For polyatomic molecules, the most relevant choice of the subsystems is not so clear cut. 
If all possible combinations (mono-, di-, tri-, . . . atomic subsystem) are included the number 
of constraints can grow quite big, and for each subsystem one needs to perform separate 
2DM optimizations. For mono-atomic subsystems this is a one-time task for a given basis 
set as all the required Hamiltonians, energies etc. can be kept stored. However, for larger 
(di-,tri-,. . . ) subsystems the problem is that the required data are geometry dependent. As 
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a consequence, for every molecular calculation one will need to also obtain the variationally 
optimized 2DM energies for all non mono-atomic subspaces. As 2DM optimizations are 
quite time consuming, this may add a lot of overhead time to a molecular 2DM optimization. 
Furthermore, the number of constraints to be included, and as such the number of extra 2DM 
energy calculations for the subspaces, grows rapidly. For a diatomic, 2 subspace constraints 
are needed, as illustrated above for BeB"''. For a molecule with 5 atoms, the number of 
constraints equals already 30 of which 25 need to be computed specifically for the molecular 
geometry considered. A further drawback is that it is sometimes hard to predict whether the 
constraint needed will involve a cationic or anionic subsystem and as a consequence, both 
are preferably included. It is very unlikely that all of these constraints are needed and only 
few will be violated when not included. Unfortunately, at this moment it is not yet clear 
how to decide a priori which are the most important constraints to be included. Despite 
these drawbacks, these constraints are very strong and their inclusion highly desirable in the 
SDP. This is shown in Ref.— where the subsystem constraints have been used in a study of 
the 14-electron isoelectronic series. 

We would like to stress that subsystem constraints are much more general then the 
atom-based constraints discussed in this paper. They hold for any subspace, e.g. also for an 
arbitrary selection of molecular orbitals (which may correspond to an active space around 
the Fermi-energy or a restricted set of orbitals of a certain symmetry). Of course, in general 
many constraints of this type will not be active. It will be interesting to see whether other 
systems can be found for which the subsystem constraints lead to improvements in the 
quality of the variationally obtained 2DM. 
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